home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / coons_warp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  7.6 KB  |  238 lines
  1. /*
  2.  * ANSI C code from the article
  3.  * "Bilinear Coons Patch Image Warping"
  4.  * by Paul S. Heckbert, ph@cs.cmu.edu
  5.  * in "Graphics Gems IV", Academic Press, 1994
  6.  *
  7.  *
  8.  * This code has been written in the most portable way possible.
  9.  * It will not compile and run without some modifications.
  10.  * You'll need to:
  11.  *    write pixel_read and pixel_write subroutines
  12.  *    read or compute a source image
  13.  *    read or compute four boundary curves
  14.  *    pick a destination image size
  15.  *    call resample() to resample the boundary curves to the width and
  16.  *        height of the destination image
  17.  *    call coons_warp() to write the destination image
  18.  *    display the warped destination image
  19.  */
  20.  
  21. #include <math.h>
  22. #include <assert.h>
  23. #define ALLOC(ptr, type, n)  assert(ptr = (type *)malloc((n)*sizeof(type)))
  24.  
  25. typedef struct {    /* 2-D POINT OR VECTOR */
  26.     float x, y;
  27. } Point2f;
  28.  
  29. typedef struct {    /* A CURVE DEFINED BY A SEQUENCE OF POINTS */
  30.     int npt;        /* number of points */
  31.     Point2f *pt;    /* array of npt points */
  32. } Curve;
  33.  
  34. #define SHIFT 20    /* number of fractional bits in fixed point coords */
  35. #define SCALE (1<<SHIFT)
  36.  
  37. typedef struct {    /* INTEGER POINT AND VECTOR */
  38.     int px, py;        /* position */
  39.     int dx, dy;        /* incremental displacement */
  40. } Ipoint;
  41.  
  42.  
  43. /*
  44.  * coons_warp: warps the picture in source image into a rectangular
  45.  * destination image according to four boundary curves, using a
  46.  * bilinear Coons patch.
  47.  * bound[0] through bound[3] are the top, right, bottom, and left
  48.  * boundary curves, respectively, clockwise from upper left.
  49.  *   (These comments are written assuming that y points down on your
  50.  *   frame buffer. Otherwise, bound should proceed CCW from lower left.)
  51.  * The lengths of bound[0] and bound[2] are assumed to be the width of
  52.  * the destination rectangle, and the lengths of bound[1] and bound[3]
  53.  * are assumed to be the height of the rectangle.
  54.  * The upper left corner of the destination rectangle is (u0,v0).
  55.  *
  56.  * Paul Heckbert    25 Feb 82, 15 Oct 93
  57.  */
  58.  
  59. void coons_warp(Pic *source, Pic *dest, Curve *bound, int u0, int v0) {
  60.     register Ipoint *pu;
  61.     register int u, x, y, qx, qy, dqx, dqy;
  62.     int nu, nv, v;
  63.     float du, dv, fv;
  64.     Point2f p00, p01, p10, p11, *pu0, *pu1, *p0v, *p1v;
  65.     Ipoint *pua;
  66.  
  67.     nu = bound[0].npt-1;            /* nu = dest_width-1 */
  68.     nv = bound[1].npt-1;            /* nv = dest_height-1 */
  69.     assert(bound[2].npt==nu+1);
  70.     assert(bound[3].npt==nv+1);
  71.  
  72.     pu0 = &bound[0].pt[0];            /* top boundary curve */
  73.     p1v = &bound[1].pt[0];            /* right */
  74.     pu1 = &bound[2].pt[nu];            /* bottom */
  75.     p0v = &bound[3].pt[nv];            /* left */
  76.     /* arrays pu1 and p0v are in the reverse of the desired order,
  77.        running from right to left and bottom to top, resp., so we
  78.        index them with negative subscripts from their ends (yeeeha!) */
  79.  
  80.     p00.x = (p0v[ 0].x + pu0[  0].x)/2.;    /* upper left patch corner */
  81.     p00.y = (p0v[ 0].y + pu0[  0].y)/2.;
  82.     p10.x = (pu0[nu].x + p1v[  0].x)/2.;    /* upper right */
  83.     p10.y = (pu0[nu].y + p1v[  0].y)/2.;
  84.     p11.x = (p1v[nv].x + pu1[-nu].x)/2.;    /* lower right */
  85.     p11.y = (p1v[nv].y + pu1[-nu].y)/2.;
  86.     p01.x = (pu1[ 0].x + p0v[-nv].x)/2.;    /* lower left */
  87.     p01.y = (pu1[ 0].y + p0v[-nv].y)/2.;
  88.  
  89.     du = 1./nu;
  90.     dv = 1./nv;
  91.  
  92.     ALLOC(pua, Ipoint, nu+1);
  93.     for (pu=pua, u=0; u<=nu; u++, pu++) {
  94.     pu->dx = (pu1[-u].x - pu0[u].x)*dv*SCALE + .5;
  95.     pu->dy = (pu1[-u].y - pu0[u].y)*dv*SCALE + .5;
  96.     pu->px = pu0[u].x*SCALE + .5;
  97.     pu->py = pu0[u].y*SCALE + .5;
  98.     }
  99.  
  100.     for (fv=0., v=0; v<=nv; v++, fv+=dv) {
  101.     qx = (p0v[-v].x - (1.-fv)*p00.x - fv*p01.x + .5)*SCALE + .5;
  102.     qy = (p0v[-v].y - (1.-fv)*p00.y - fv*p01.y + .5)*SCALE + .5;
  103.     dqx = (p1v[v].x - p0v[-v].x - (1.-fv)*(p10.x-p00.x) - fv*(p11.x-p01.x))
  104.         *du*SCALE + .5;
  105.     dqy = (p1v[v].y - p0v[-v].y - (1.-fv)*(p10.y-p00.y) - fv*(p11.y-p01.y))
  106.         *du*SCALE + .5;
  107.     for (pu=pua, u=0; u<=nu; u++, pu++) {
  108.         x = pu->px+qx >> SHIFT;
  109.         y = pu->py+qy >> SHIFT;
  110.         pixel_write(dest, u0+u, v0+v, pixel_read(source, x, y));
  111.         qx += dqx;
  112.         qy += dqy;
  113.         pu->px += pu->dx;
  114.         pu->py += pu->dy;
  115.     }
  116.     }
  117.     free(pua);
  118. }
  119.  
  120. /* the following routine is the slow way to do a Coons warp, for reference */
  121.  
  122. void coons_warp1(Pic *source, Pic *dest, Curve *bound, int u0, int v0) {
  123.     int nu, nv, u, v, x, y;
  124.     float fu, fv;
  125.     Point2f p00, p01, p10, p11, *pu0, *pu1, *p0v, *p1v;
  126.  
  127.     nu = bound[0].npt-1;
  128.     nv = bound[1].npt-1;
  129.  
  130.     pu0 = &bound[0].pt[0];            /* top boundary curve */
  131.     p1v = &bound[1].pt[0];            /* right */
  132.     pu1 = &bound[2].pt[nu];            /* bottom */
  133.     p0v = &bound[3].pt[nv];            /* left */
  134.  
  135.     p00.x = (p0v[ 0].x + pu0[  0].x)/2.;    /* upper left patch corner */
  136.     p00.y = (p0v[ 0].y + pu0[  0].y)/2.;
  137.     p10.x = (pu0[nu].x + p1v[  0].x)/2.;    /* upper right */
  138.     p10.y = (pu0[nu].y + p1v[  0].y)/2.;
  139.     p11.x = (p1v[nv].x + pu1[-nu].x)/2.;    /* lower right */
  140.     p11.y = (p1v[nv].y + pu1[-nu].y)/2.;
  141.     p01.x = (pu1[ 0].x + p0v[-nv].x)/2.;    /* lower left */
  142.     p01.y = (pu1[ 0].y + p0v[-nv].y)/2.;
  143.  
  144.     for (v=0; v<=nv; v++) {
  145.     fv = (float)v/nv;
  146.     for (u=0; u<=nu; u++) {
  147.         fu = (float)u/nu;
  148.         x =     (1.-fv)*pu0[ u].x + fv*pu1[-u].x
  149.            + (1.-fu)*p0v[-v].x + fu*p1v[ v].x
  150.            - (1.-fu)*(1.-fv)*p00.x - fu*(1.-fv)*p10.x
  151.            - (1.-fu)*   fv    *p01.x - fu*   fv  *p11.x + .5;
  152.         y =     (1.-fv)*pu0[ u].y + fv*pu1[-u].y
  153.            + (1.-fu)*p0v[-v].y + fu*p1v[ v].y
  154.            - (1.-fu)*(1.-fv)*p00.y - fu*(1.-fv)*p10.y
  155.            - (1.-fu)*   fv    *p01.y - fu*   fv  *p11.y + .5;
  156.         pixel_write(dest, u0+u, v0+v, pixel_read(source, x, y));
  157.     }
  158.     }
  159. }
  160.  
  161. /* resample_nonuniform: non-uniformly resample curve a to create curve b,
  162.  * with n points.  Allocates b->pt to have length n */
  163.  
  164. void resample_nonuniform(Curve *a, Curve *b, int n) {
  165.     int ai, bi;
  166.     double ax, af;
  167.     Point2f *ap, *bp;
  168.  
  169.     if (n<2) {
  170.     fprintf(stderr, "Only %d point in new curve\n", n);
  171.     exit(1);
  172.     }
  173.     ALLOC(b->pt, Point2f, n);
  174.     b->npt = n;
  175.     for (bp=b->pt, bi=0; bi<n; bi++, bp++) {
  176.     ax = (float)bi*(a->npt-1)/(n-1);
  177.     ai = ax;
  178.     af = ax-ai;
  179.     ap = &a->pt[ai];
  180.     bp->x = ap[0].x + af*(ap[1].x-ap[0].x);
  181.     bp->y = ap[0].y + af*(ap[1].y-ap[0].y);
  182.     }
  183. }
  184.  
  185. static double len(double x, double y) {return sqrt(x*x+y*y);}
  186.  
  187. /* resample_uniform: uniformly resample curve a to create curve b,
  188.  * with n points.  Allocates b->pt to have length n */
  189.  
  190. void resample_uniform(Curve *a, Curve *b, int n) {
  191.     int i;
  192.     double step, l, d;
  193.     Point2f *ap, *bp;
  194.  
  195.     if (a->npt<2) {
  196.     fprintf(stderr, "Only %d point in curve\n", a->npt);
  197.     exit(1);
  198.     }
  199.     for (step=0., ap=a->pt, i=a->npt-1; i>0; i--, ap++)
  200.     step += len(ap[1].x-ap[0].x, ap[1].y-ap[0].y);
  201.     step /= n-1;        /* length of each output segment (ideally) */
  202.     ALLOC(b->pt, Point2f, n);
  203.     b->npt = n;
  204.     d = .0001;            /* = 0 + tolerance for roundoff error */
  205.     for (ap=a->pt, bp=b->pt, i=a->npt-1; i>0; i--, ap++) {
  206.     l = len(ap[1].x-ap[0].x, ap[1].y-ap[0].y);
  207.     d += l;
  208.     /* d is the remaining length of the line segment from ap[0] to ap[1]
  209.      * that needs to be subdivided into segments of length step */
  210.     while (d>0.) {
  211.         bp->x = ap[1].x - d/l*(ap[1].x-ap[0].x);
  212.         bp->y = ap[1].y - d/l*(ap[1].y-ap[0].y);
  213.         bp++;
  214.         d -= step;
  215.     }
  216.     }
  217.     if (bp-b->pt != n)
  218.     printf("WARNING: requested %d points, created %d, d=%g\n",
  219.         n, bp-b->pt, d);
  220. }
  221.  
  222. static void resample(Curve *a, Curve *b, int n, int nonuniform) {
  223.     if (nonuniform) resample_nonuniform(a, b, n);
  224.     else resample_uniform(a, b, n);
  225. }
  226.  
  227. Curve *boundary_resample(Curve *a, int nx, int ny, int nonuniform) {
  228.     /* "nonuniform" is a boolean flag */
  229.     Curve *b;
  230.  
  231.     ALLOC(b, Curve, 4);
  232.     resample(&a[0], &b[0], nx, nonuniform);
  233.     resample(&a[1], &b[1], ny, nonuniform);
  234.     resample(&a[2], &b[2], nx, nonuniform);
  235.     resample(&a[3], &b[3], ny, nonuniform);
  236.     return b;
  237. }
  238.